{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Kittens"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Modeling and Simulation in Python*\n",
    "\n",
    "Copyright 2021 Allen Downey\n",
    "\n",
    "License: [Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International](https://creativecommons.org/licenses/by-nc-sa/4.0/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# download modsim.py if necessary\n",
    "\n",
    "from os.path import basename, exists\n",
    "\n",
    "def download(url):\n",
    "    filename = basename(url)\n",
    "    if not exists(filename):\n",
    "        from urllib.request import urlretrieve\n",
    "        local, _ = urlretrieve(url, filename)\n",
    "        print('Downloaded ' + local)\n",
    "    \n",
    "download('https://github.com/AllenDowney/ModSimPy/raw/master/modsim.py')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# import functions from modsim\n",
    "\n",
    "from modsim import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you have used the Internet, you have probably seen videos of kittens unrolling toilet paper.\n",
    "And you might have wondered how long it would take a standard kitten to unroll 47 m of paper, the length of a standard roll.\n",
    "\n",
    "The interactions of the kitten and the paper rolls are complex.  To keep things simple, let's assume that the kitten pulls down on the free end of the roll with constant force.  And let's neglect the friction between the roll and the axle.\n",
    "\n",
    "This diagram shows the paper roll with the force applied by the kitten, $F$, the lever arm of the force around the axis of rotation, $r$, and the resulting torque, $\\tau$.\n",
    "\n",
    "![Diagram of a roll of toilet paper, showing a force, lever arm, and the resulting torque.](https://github.com/AllenDowney/ModSim/raw/main/figs/kitten.png)\n",
    "\n",
    "Assuming that the force applied by the kitten is 0.002 N, how long would it take to unroll a standard roll of toilet paper?\n",
    "\n",
    "We'll use the same parameters as in Chapter 24:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "Rmin = 0.02      # m\n",
    "Rmax = 0.055     # m\n",
    "Mcore = 15e-3    # kg\n",
    "Mroll = 215e-3   # kg\n",
    "L = 47           # m\n",
    "tension = 0.002  # N"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Rmin` and `Rmax` are the minimum and maximum radius of the roll, respectively.\n",
    "`Mcore` is the weight of the core (the cardboard tube at the center) and `Mroll` is the total weight of the paper.\n",
    "`L` is the unrolled length of the paper.\n",
    "`tension` is the force the kitten applies by pulling on the loose end of the roll (I chose this value because it yields reasonable results)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In Chapter 24 we defined $k$ to be the constant that relates a change in the radius of the roll to a change in the rotation of the roll:\n",
    "\n",
    "$$dr = k~d\\theta$$ \n",
    "\n",
    "And we derived the equation for $k$ in terms of $R_{min}$, $R_{max}$, and $L$. \n",
    "\n",
    "$$k =  \\frac{1}{2L} (R_{max}^2 - R_{min}^2)$$\n",
    "\n",
    "So we can compute `k` like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = (Rmax**2 - Rmin**2) / 2 / L    \n",
    "k    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Moment of Inertia\n",
    "\n",
    "To compute angular acceleration, we'll need the moment of inertia for the roll.\n",
    "\n",
    "At <http://modsimpy.com/moment> you can find moments of inertia for\n",
    "simple geometric shapes. I'll model the core as a \"thin cylindrical shell\", and the paper roll as a \"thick-walled cylindrical tube with open ends\".\n",
    "\n",
    "The moment of inertia for a thin shell is just $m r^2$, where $m$ is the mass and $r$ is the radius of the shell.\n",
    "\n",
    "For a thick-walled tube the moment of inertia is\n",
    "\n",
    "$$I = \\frac{\\pi \\rho h}{2} (r_2^4 - r_1^4)$$ \n",
    "\n",
    "where $\\rho$ is the density of the material, $h$ is the height of the tube (if we think of the roll oriented vertically), $r_2$ is the outer diameter, and $r_1$ is the inner diameter.\n",
    "\n",
    "Since the outer diameter changes as the kitten unrolls the paper, we\n",
    "have to compute the moment of inertia, at each point in time, as a\n",
    "function of the current radius, `r`, like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def moment_of_inertia(r):\n",
    "    \"\"\"Moment of inertia for a roll of toilet paper.\n",
    "    \n",
    "    r: current radius of roll in meters\n",
    "    \n",
    "    returns: moment of inertia in kg m**2\n",
    "    \"\"\"    \n",
    "    Icore = Mcore * Rmin**2   \n",
    "    Iroll = np.pi * rho_h / 2 * (r**4 - Rmin**4)\n",
    "    return Icore + Iroll"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Icore` is the moment of inertia of the core; `Iroll` is the moment of inertia of the paper.\n",
    "\n",
    "`rho_h` is the density of the paper in terms of mass per unit of area. \n",
    "To compute `rho_h`, we compute the area of the complete roll like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "area = np.pi * (Rmax**2 - Rmin**2)\n",
    "area"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And divide the mass of the roll by that area."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "rho_h = Mroll / area\n",
    "rho_h"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As an example, here's the moment of inertia for the complete roll."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "moment_of_inertia(Rmax)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As `r` decreases, so does `I`.  Here's the moment of inertia when the roll is empty."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "moment_of_inertia(Rmin)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The way $I$ changes over time might be more of a problem than I have made it seem.  In the same way that $F = m a$ only applies when $m$ is constant, $\\tau = I \\alpha$ only applies when $I$ is constant.  When $I$ varies, we usually have to use a more general version of Newton's law.  However, I believe that in this example, mass and moment of inertia vary together in a way that makes the simple approach work out.\n",
    "\n",
    "A friend of mine who is a physicist is not convinced; nevertheless, let's proceed on the assumption that I am right."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulation\n",
    "\n",
    "The state variables we'll use are\n",
    "\n",
    "* `theta`, the total rotation of the roll in radians, \n",
    "\n",
    "* `omega`, angular velocity in rad / s,\n",
    "\n",
    "* `r`, the radius of the roll, and\n",
    "\n",
    "* `y`, the length of the unrolled paper.\n",
    "\n",
    "Here's a `State` object with the initial conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "init = State(theta=0, omega=0, y=0, r=Rmax)\n",
    "init"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here's a `System` object with the starting conditions and `t_end`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "system = System(init=init, t_end=120)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can take it from here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:**\n",
    "\n",
    "Write a slope function we can use to simulate this system.  Test it with the initial conditions.  The results should be approximately\n",
    "\n",
    "```\n",
    "0.0, 0.294, 0.0, 0.0\n",
    "```\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Write an event function that stops the simulation when `y` equals `L`, that is, when the entire roll is unrolled.  Test your function with the initial conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now run the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And check the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "results.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The final value of `theta` should be about 200 rotations, the same as in Chapter 24.\n",
    "\n",
    "The final value of `omega` should be about 63 rad/s, which is about 10 revolutions per second.  That's pretty fast, but it might be plausible.\n",
    "\n",
    "The final value of `y` should be `L`, which is 47 m.\n",
    "\n",
    "The final value of `r` should be `Rmin`, which is 0.02 m.\n",
    "\n",
    "And the total unrolling time should be about 76 seconds, which seems plausible."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cells plot the results.\n",
    "\n",
    "`theta` increases slowly at first, then accelerates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "results.theta.plot(color='C0', label='theta')\n",
    "decorate(xlabel='Time (s)',\n",
    "         ylabel='Angle (rad)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Angular velocity, `omega`, increases almost linearly at first, as constant force yields almost constant torque.  Then, as the radius decreases, the lever arm decreases, yielding lower torque, but moment of inertia decreases even more, yielding higher angular acceleration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "results.omega.plot(color='C2', label='omega')\n",
    "\n",
    "decorate(xlabel='Time (s)',\n",
    "         ylabel='Angular velocity (rad/s)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`y` increases slowly and then accelerates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "results.y.plot(color='C1', label='y')\n",
    "\n",
    "decorate(xlabel='Time (s)',\n",
    "         ylabel='Length (m)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`r` decreases slowly, then accelerates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "results.r.plot(color='C4', label='r')\n",
    "\n",
    "decorate(xlabel='Time (s)',\n",
    "         ylabel='Radius (m)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
